function q = output_corner(mu,B,A,P,rho,lambda_rnd,lambda_tech,s)

global options_fmincon;
global options_QP;
global options_SA;

opt_alg = 'qp'; % 'qp' 'sa' 'fmincon'

%% Compute trial analytic interior solution.
n = length(A);
nu = mu + (eye(n) + lambda_rnd*A + lambda_tech*P)*s;
H = eye(n) + rho*B - lambda_rnd*A - lambda_tech*P; % Hessian.
f = -nu; 
M = inv(H);
R = M*(eye(n) + lambda_rnd*A + lambda_tech*P);
q_bar = M*mu;
q = M*mu + R*s;

%% If trial analytic interior solution has negative output levels, compute trial analytic corner solution, by setting all firms with negative output levels in the trial analytic interior solution to zero.
ind = find(q<0);
if(~isempty(ind))

    ind_corner = find(q<=0);
    ind_interior = find(q>0);
    
    A_interior_interior = A(ind_interior,ind_interior); 
    A_interior_corner = A(ind_interior,ind_corner);
    A_corner_interior = A(ind_corner,ind_interior);
    A_corner_corner = A(ind_corner,ind_corner);
    B_corner_interior = B(ind_corner,ind_interior);
    P_interior_interior = P(ind_interior,ind_interior);
    P_interior_corner = P(ind_interior,ind_corner);
    P_corner_interior = P(ind_corner,ind_interior);
    P_corner_corner = P(ind_corner,ind_corner);
    M_interior_interior = M(ind_interior,ind_interior);
    R_interior_interior = R(ind_interior,ind_interior);
    s_interior = s(ind_interior);
    s_corner = s(ind_corner);
    mu_interior = mu(ind_interior);
    mu_corner = mu(ind_corner);
    
    q_sub_interior = M_interior_interior*mu_interior + R_interior_interior*s_interior ...
        + lambda_rnd*A_interior_corner*s_corner + lambda_tech*P_interior_corner*s_corner;
    
    cond_corner = - mu_corner - s_corner ...
        + rho*B_corner_interior*q_sub_interior ...
        - lambda_rnd*A_corner_interior*s_interior - lambda_tech*P_corner_interior*s_interior ...
        - lambda_rnd*A_corner_corner*s_corner - lambda_tech*P_corner_corner*s_corner;
    
    q = zeros(n,1);
    for i=1:length(ind_interior)
        q(ind_interior(i)) = q_sub_interior(i);
    end
    
    %% If trial analytic corner solution violates corner solution conditions, compute corner solution using an optimization algorithm with output levels
    %% from above as starting values.
    ind = find(q_sub_interior<0);
    ind2 = find(cond_corner<0);
    if(~isempty(ind) || ~isempty(ind2))

        lb = zeros(n,1);
        fun = @(x) (1/2*(x'*H*x) + f*x);
        switch opt_alg
            case 'fmincon'
                q = fmincon(fun,q,[],[],[],[],lb,[],[],options_fmincon);
            case 'qp'
                q = quadprog(H,f,[],[],[],[],lb,[],q,options_QP);
            case 'sa'
                [q,fval,exitFlag,out] = simulannealbnd(fun,q,lb,[],options_SA);
        end

        ind_corner = find(q<=0);
        ind_interior = find(q>0);
        
        A_interior_interior = A(ind_interior,ind_interior);
        A_interior_corner = A(ind_interior,ind_corner);
        A_corner_interior = A(ind_corner,ind_interior);
        A_corner_corner = A(ind_corner,ind_corner);
        B_corner_interior = B(ind_corner,ind_interior);
        P_interior_interior = P(ind_interior,ind_interior);
        P_interior_corner = P(ind_interior,ind_corner);
        P_corner_interior = P(ind_corner,ind_interior);
        P_corner_corner = P(ind_corner,ind_corner);
        M_interior_interior = M(ind_interior,ind_interior);
        R_interior_interior = R(ind_interior,ind_interior);
        q_sub_interior = q(ind_interior);
        s_interior = s(ind_interior);
        s_corner = s(ind_corner);
        mu_interior = mu(ind_interior);
        mu_corner = mu(ind_corner);
        
        cond_corner = - mu_corner - s_corner ...
            + rho*B_corner_interior*q_sub_interior ...
            - lambda_rnd*A_corner_interior*s_interior - lambda_tech*P_corner_interior*s_interior ...
            - lambda_rnd*A_corner_corner*s_corner - lambda_tech*P_corner_corner*s_corner;
        
        ind = find(cond_corner<0);
        if(~isempty(ind2))
            disp(['Corner solution condition violated for numerically computed Nash output levels!'])
        end
        
    end
    
end

%% Check if the Hessian is positive definite.
if(~isempty(find(eig(H)<0)))
    
    disp(['Hessian is not positive definite. Search for multiple equilibria using simulated annealing algorithm...'])
    
    %% Run simulated annealing algorithm repeatedly.
    lb = zeros(n,1);
    fun = @(x) (1/2*(x'*H*x) + f*x);
    for k=1:20
        [q_k,fval_k,exitFlag,out] = simulannealbnd(fun,q,lb,[],options_SA);
        q = [q: q_k];
    end
    
end

